clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"

cap log close
log using $figures\J_file,replace t

cd $data

************************************************************************************************************************************************************************************
******VIGNETTE DATA
************************************************************************************************************************************************************************************
************************************************************************************************************************************************************************************
u $data\vignette_data, clear
******************************************************************************************************************************************************************

gen 	vign_sev_f=vign_severe 
replace vign_sev_f=. if vign_f==0
gen 	vign_sev_m=vign_severe
replace vign_sev_m=. if vign_f==1

ren walkra walkra_R_1 
ren dressa dressa_R_1 
ren stoopa stoopa_R_1 
ren beda   beda_R_1   
ren bmi    bmi_R      
ren hosp   hosp_R      

qui tab occ_interm,gen(occx)
qui tab vign_id,gen(vigndomx)

label var vign_severe "Overall"
label var vign_sev_f "female hypoth. person"
label var vign_sev_m "male hypoth. person"
label var cantwork "Respondent disabled"
label var age "Respondent's age"

lab var doctor_told_has_hbp "Have HBP" 
lab var doctor_told_has_psy "Have psych. cond."
lab var doctor_told_has_hea "Have heart cond."
lab var doctor_told_has_art "Have arthritis"
lab var doctor_told_has_dia "Have diabetes"
lab var doctor_told_has_lun "Have lung dis."
lab var doctor_told_has_str "Have stroke"
lab var doctor_told_has_can "Have cancer"
label var vign_severe "Overall"
label var vign_sev_f "female hypoth. person"
label var vign_sev_m "male hypoth. person"
label var cantwork "Respondent disabled"
label var age "Respondent's age"
label var female "Female respondent"
label var vign_f "Female vignette"
label var age "Age"
label var college "College degree"
label var racex2 "Black"
lab var walkra "ADL: walking"
lab var dressa "ADL: dressing"
lab var stoopa "ADL: stooping"
lab var beda "ADL: getting out of bed" 
lab var hosp "Spent any night in hosp" 
lab var bmi  "BMI" 
lab def occ_interm 1 "Unknown" 2 "Manag/Prof" 3 "Sales/clerical" 4 "Clean/Protect/Food serv" 5 "Personal/Health serv" ///
	6 "Farm/Constr/Mech" 7 "Precision/Armed force" 8 "Operators"
lab values occ_interm occ_interm
lab def vign_id 1 "Low pain" 2 "Interm. pain" 3 "Severe pain" ///
				4 "Low depression" 5 "Interm. depression" 6 "Severe depression" ///
				7 "Low cardiov" 8 "Interm. cardiov" 9 "Severe cardiov"
lab values vign_id vign_id
lab var racex2 "Black"
lab var married "Married" 
lab var widowed "Widowed"


lab var occx1 "Unknown" 
lab var occx2 "Manag/Prof" 
lab var occx3 "Sales/clerical" 
lab var occx4 "Clean/Protect/Food serv" 
lab var occx5 "Personal/Health serv" 
lab var occx6 "Farm/Constr/Mech"
lab var occx7 "Precision/Armed force" 
lab var occx8 "Operators"

label var female "Female"
label var age "Age"
label var college "College degree"
label var racex2 "Black"

************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_L 		female college racex2 doctor* married widowed age  walkra dressa stoopa beda hosp bmi i.occ_interm experience
global spec_V 		$spec_L vign_f i.vign_id 
************************************************************************************************************************************************************************************

******************************************************************************************************************************************************************
**Sample selection (drop those w/ missing data in the variables used for the 2 regressions)
qui biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	
keep if e(sample)
******************************************************************************************************************************************************************

sort hhidpn vign_id
gen cw=cantwork 
replace cw=. if hhidpn==hhidpn[_n-1]
gen ag=age 
replace ag=. if hhidpn==hhidpn[_n-1]

lab var cw "Respondent disabled"
lab var ag "Respondent's age"

eststo: estpost sum vign_severe vign_sev_m vign_sev_f cw ag
bysort female: eststo: estpost sum vign_severe vign_sev_m vign_sev_f cw ag

esttab using $figures\tableOA3.tex, cells("mean(fmt(%9.2f)) sd(fmt(%9.2f))") label nodepvar replace ///
 title("Vignettes: Descriptive Statistics") ///
mtitle("All resp." "Male resp." "Female resp.")
eststo clear

drop cw ag

label var female "Female respondent"
label var vign_f "Female vignette"

hetprobit vign_severe $spec_V ,het(female)
scalar het_V=_b[lnsigma:female]
scalar het_se_V=_se[lnsigma:female]

	
*This is the bivariate probit L(respondent), L(vignette)
****TABLE 9 + \RHO_{L,V}
eststo bigml2: biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	

matrix bV=e(b)
matrix bV=bV'

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments	=_b[vign_severe:female]\_b[vign_severe:vign_f]\r(table)[1,colsof(r(table))]
matrix semoments=_se[vign_severe:female]\_se[vign_severe:vign_f]\r(table)[2,colsof(r(table))]	
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
************************************************************************************************************************************************

************************************************************************************************************************************************
estadd scalar corr_L_V 	r(table)[1,colsof(r(table))]

qui su vign_severe
scalar sample_avg_V=r(mean)
scalar beta_mom_V=_b[vign_severe:female]
scalar beta_mom2_V=_b[vign_severe:vign_f]

estadd scalar het_V
estadd scalar het_se_V

************************************************************************************************************************************************


************************************************************************************************************************************************************************************
********REJECTION DATA *************************************************************************************************************************************************************************************
*************************************************************************************************************************************************************************************

u $data\rejection_data, clear

/*This is -X*a(Lbar)*/
#delimit ;
gen a1_fromV=
bV[2 ,1]		*college+  
bV[3 ,1]         *racex2  +
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[12,1]         *married +
bV[13,1]         *widowed +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +
bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr

/*This is -X*a(Lbar)(institutional variables only)*/
#delimit ;
gen a1_inst_fromV=
bV[2 ,1]		*college+  
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +
bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr

************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global c_L 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm 
global c_A 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm logbs 
global c_R 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9
************************************************************************************************************************************************************************************
gen A=applied
gen R=nosuccess
gen L=cantwork_notemp

label var L "Work lim."
label var age "Age"
lab var bsx1 "BS=Musculoskeletal"
lab var bsx2 "BS=Respiratory"   
lab var bsx3 "BS=Cardiov."   
lab var bsx4 "BS=Endocrine"       
lab var bsx5 "BS=Neurol."
lab var bsx6 "BS=Mental dis."
lab var bsx7 "BS=Cancer"    
lab var bsx8 "BS=Immune def."  
lab var bsx9 "BS=Dig. and Urin." 
lab var bsx10 "BS=Other"
lab var doctor_told_has_hbp "Have HBP" 
lab var doctor_told_has_psy "Have psych. cond."
lab var doctor_told_has_hea "Have heart cond."
lab var doctor_told_has_art "Have arthritis"
lab var doctor_told_has_dia "Have diabetes"
lab var doctor_told_has_lun "Have lung dis."
lab var doctor_told_has_str "Have stroke"
lab var doctor_told_has_can "Have cancer"
label var experience "Experience"
label var female "Female"
label var college "College degree"
label var racex2 "Black"
lab var walkra "ADL: walking"
lab var dressa "ADL: dressing"
lab var stoopa "ADL: stooping"
lab var beda "ADL: getting out of bed" 
lab var hosp "Spent any night in hosp" 
lab var bmi  "BMI" 
lab def occ_interm 1 "Unknown" 2 "Manag/Prof" 3 "Sales/cler." 4 "Clean/Protect/Food serv" 5 "Pers./Health serv" ///
	6 "Farm/Constr/Mech" 7 "Precision/Army" 8 "Operators",replace
lab values occ_interm occ_interm
lab var married "Married" 
lab var widowed "Widowed"
lab var A "Applied"
lab var R "Reject."
lab var logbs "Log(liquid res.)"

****HETEROSKEDASTICITY TEST for L
qui hetprobit L $c_L ,het(female)
scalar 			 het_L=_b[lnsigma:female]
scalar 		  het_se_L=_se[lnsigma:female]
****HETEROSKEDASTICITY TEST for A
qui hetprobit A $c_A ,het(female)
scalar 			 het_A=_b[lnsigma:female]
scalar 		  het_se_A=_se[lnsigma:female]
****HETEROSKEDASTICITY TEST for R|A
cap program drop heckprobithet
program heckprobithet
	version 17.0
	args lnf theta1 theta2 theta3
	tempvar p
	qui gen double `p'=(1+exp(2*`theta3'))^(-0.5)
	quietly replace `lnf'=ln(binormal(`p'*`theta1',`theta2',`p')) 	if $ML_y1==1
	quietly replace `lnf'=ln(binormal(-`p'*`theta1',`theta2',-`p')) if $ML_y1==2
	quietly replace `lnf'=ln(1-normal(`theta2')) 					if $ML_y1==3
end

gen 	index2=1 if R==1 & A==1
replace index2=2 if R==0 & A==1
replace index2=3 if A==0
ml model lf heckprobithet (index2 = $c_R) (theta2: $c_A) (theta3: female,nocons)
qui ml maximize, nolog
scalar het_R = _b[theta3:female]
scalar het_se_R=_se[theta3:female]

set seed 2134569

// Create Halton draw variables
mdraws, dr(10) neq(3) prefix(z)			/*If dr(N), then in egen `sp3' and egen `sp2' adjust dr(N) accordingly*/

// Get initial values
quietly {
        probit A $c_A
        mat b1 = e(b)
        mat coleq b1 = A
        probit L $c_L
        mat b2 = e(b)
        mat coleq b2 = L
        probit R $c_R
        mat b3 = e(b)
        mat coleq b3 = R
        mat b0 = b1, b2, b3
}


cap program drop myll2
program define myll2
args lnf xb1 xb2 xb3 c21 c31 c32
        tempvar sp2 sp3 k1 k2 k3
        quietly {
            gen double `k1' = 2*$ML_y1 - 1
            gen double `k2' = 2*$ML_y2 - 1
            gen double `k3' = 2*$ML_y3 - 1
            tempname cf21 cf31 cf32 cf22 cf33 C1 C2
            su `c21', meanonly
            scalar `cf21' = r(mean)
            su `c31', meanonly
            scalar `cf31' = r(mean)
            su `c32', meanonly
            scalar `cf32' = r(mean)
            scalar `cf22' = sqrt( 1 - `cf21'^2 )
            scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 )
            tempname C1 C2
            mat `C1' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31', `cf32', `cf33')
            mat `C2' = (1, 0 \ `cf21', `cf22' )
            egen `sp3' = mvnp(`xb1' `xb2' `xb3') if $ML_y1==1, ///
                    chol(`C1') dr(10) prefix(z) signs(`k1' `k2' `k3') adoonly
            egen `sp2' = mvnp(`xb1' `xb2' ) if $ML_y1==0, ///
                    chol(`C2') dr(10) prefix(z) signs(`k1' `k2') adoonly
            replace `lnf'= ln(`sp3') if $ML_y1==1
            replace `lnf'= ln(`sp2') if $ML_y1==0
        }
end

ml model lf myll2 	(A: A=$c_A) ///
					(L: L=$c_L) ///
					(R: R=$c_R) ///
					/c21 /c31 /c32 ///
					, missing

ml init b0
eststo bigml1:ml maximize
qui su A
scalar sample_avg_A=r(mean)
qui su L
scalar sample_avg_L=r(mean)
qui su R
scalar sample_avg_R=r(mean)

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[L:female]\_b[A:female]
matrix semoments=semoments\_se[L:female]\_se[A:female]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar betaFL=_b[L:female]	
scalar delta =_b[R:female]	/*to be used for pinning down \alpha*/

gen bF_L=_b[L:female]
gen bF_A=_b[A:female] 
gen bF_R=_b[R:female]
	
#delimit;
gen a1_age=_b[L:age]*age;
gen a1_institutional_year=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R+
_b[L:age]*                           age+
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8)+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15);
gen a1_health=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R;
gen a1_demo=
_b[L:female]*				      female+
_b[L:age]*                           age+
_b[L:racex2]*                     racex2+
_b[L:married]*                   married+
_b[L:widowed]*                   widowed+
_b[L:_cons];
gen a1_labor=
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15)+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a1=a1_demo+a1_health+a1_labor
foreach x in a1 a1_institutional_year a1_demo a1_health a1_labor a1_age	{
	replace `x'=. if L==.
}

#delimit;
gen a2_age=_b[A:age]*age;
gen a2_institutional_year=
_b[A:age]*                           age+
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R+
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8)+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15);
gen a2_demo=
_b[A:female]*				      female+
_b[A:age]*                           age+
_b[A:racex2]*                     racex2+
_b[A:married]*                   married+
_b[A:widowed]*                   widowed+
_b[A:logbs]*logbs						+
_b[A:_cons];
gen a2_health=
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R;
gen a2_labor=
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15)+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a2=a2_demo+a2_health+a2_labor
foreach x in a2 a2_institutional_year a2_demo a2_health a2_labor a2_age	{
	replace `x'=. if A==.
}

#delimit;
gen a3_age=_b[R:age]*age;
gen a3_institutional_year=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9+
	_b[R:age]*                           age+
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8)+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020);	
gen a3_health=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9;
gen a3_demo=
	_b[R:female]*				      female+
	_b[R:age]*                           age+
	_b[R:racex2]*                     racex2+
	_b[R:married]*                   married+
	_b[R:widowed]*                   widowed+
	_b[R:_cons];
gen a3_labor=
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020)+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen a3=a3_demo+a3_health+a3_labor

foreach x in a3 a3_institutional_year a3_demo a3_health a3_labor a3_age	{
	replace `x'=. if R==.
}

nlcom (r21: /c21) (r31: /c31) (r32: /c21*/c31 + sqrt( 1 - /c21^2 )*/c32)
estadd scalar corr_A_L r(b)[1,1]
estadd scalar corr_A_R r(b)[1,2]
estadd scalar corr_L_R r(b)[1,3]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\r(b)[1,1]\r(b)[1,2]\r(b)[1,3]
matrix semoments=semoments\sqrt(r(V)[1,1])\sqrt(r(V)[2,2])\sqrt(r(V)[3,3])
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_aw_ap=-r(b)[1,2]	/*It's (A,L,R), so r(1,1)=corr(A,L) and r(1,2)=corr(A,R), so corr(Aw,A)=-r(1,2) */

gen rho_A_L=r(b)[1,1]
gen rho_A_R=r(b)[1,2]
gen rho_L_R=r(b)[1,3]

preserve
	/*The vignette identifies -X*a(Lbar)=a1_fromV	*/
	/*The work limitations identifies X*a(L*)-X*aL(bar)=a1 (from which need to subtract the female effect) */
	/*To retrieve X*a(L*)=a1+X*aL(bar)=a1-a1_fromV. That's what we have below*/
	gen xbL=(a1-betaFL*female)-(a1_fromV)			 
	gen xbL_only_inst=(a1-betaFL*female)-(a1_inst_fromV)			 
	keep hhidpn a1 a2 xbL* female A R college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R year occ_interm bsx1-bsx9 
	save $data\data_to_est_alpha,replace			
	/* Do the saving here, before dropping the first 2 waves for estimating log spending */
restore

qui tab wave,gen(waved) 
qui tab occ_interm,gen(occx)
qui tab year,gen(yeard)

save $data\after_ML,replace




u $data\after_ML,clear
drop rho* bF*

*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_C 	female college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave ///
				i.occ_interm logy lowcash govt_hins priv_hins 		
global spec_U 	$spec_C spin somecov* sp_medexp*		
************************************************************************************************************************************************************************************

gen Aw=1-R

gen A_Aw_3ml=-a3
gen A_App_3ml=a2


*********************************************************
drop if wave<3 			/*no spending data in waves 1,2*/
*********************************************************

duplicates tag hhidpn year logx,gen(check_type12_only)			/*Need 1 obs per year/per HH*/
tab check_type12_only,miss

gen nonzero=x>0
gen uncensored_C=logx!=.

***TABLE 10 + SIGMA_{C}  
eststo bigml3: heckman logx $spec_C, select(nonzero= $spec_U ) nolog
gen logxsample=e(sample)
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[logx:female]\r(table)[1,colsof(r(table))-1]
matrix semoments=semoments\_se[logx:female]\r(table)[2,colsof(r(table))-1]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
scalar sigmac	=r(table)[1,colsof(r(table))-1]
scalar se_sigmac=r(table)[2,colsof(r(table))-1]
qui su logx
scalar sample_avg_x=r(mean)

predict xb,xb
gen u=logx-xb

lab var lowcash	"Less than USD 500"
lab var logy "Log(income)"
lab var govt_hins "Govt health ins"
lab var priv_hins "Priv health ins"
lab var sp_medexp   "Spouse with pos. OOP" 
lab var sp_medexp2  "Spouse's OOP" 
lab var spin "Spouse has h. ins"
lab var somecov1 "Full covg of docs,etc"
lab var somecov2 "Full covg of hosp,etc"
estadd scalar het_L ,replace
estadd scalar het_se_L ,replace
estadd scalar het_A ,replace
estadd scalar het_se_A ,replace
estadd scalar het_R ,replace
estadd scalar het_se_R,replace
estadd scalar sigmac,replace
estadd scalar se_sigmac,replace

estadd scalar sample_avg_V, replace
estadd scalar sample_avg_A, replace
estadd scalar sample_avg_L, replace
estadd scalar sample_avg_R, replace
estadd scalar sample_avg_x, replace

drop waved* occx* yeard*

biprobit (uncensored_C =$spec_U) (L =	$c_L),nolog	
predict a_unc1,xb1
predict a_L,xb2
scalar rho_uC_uL=e(rho)

gen corr_C_0=normalden(a_unc1)*(normal((1-rho_uC_uL^2)^(-0.5)*(a_L-rho_uC_uL*a_unc1)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)		/*These expressions are the ones verified via Monte Carlo simulations*/
gen corr_C_L=normalden(a_L)   *(normal((1-rho_uC_uL^2)^(-0.5)*(a_unc1-rho_uC_uL*a_L)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)

***TABLE 10: RHO_{C,L}
reg u corr_C_0 corr_C_L if L==1 & uncensored_C==1,nocons
gen rcl = _b[corr_C_L]
 
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_L]						/*This is corr(u_C,u_L)*/
matrix semoments=semoments\_se[corr_C_L]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
											
												
biprobit (uncensored_C =$spec_U) (A =	$c_A),nolog	
predict a_unc2,xb1
predict a_App,xb2
scalar rho_uC_uA=e(rho)
scalar rho_unc_ap=e(rho)
	
gen corr_C_0_2=normalden(a_unc2)*(normal((1-rho_uC_uA^2)^(-0.5)*(a_App-rho_uC_uA*a_unc2)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)	
/*These expressions have been verified via MC simulations*/
gen corr_C_App=normalden(a_App) *(normal((1-rho_uC_uA^2)^(-0.5)*(a_unc2-rho_uC_uA*a_App)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)

***TABLE 10: RHO_{C,A}
reg u corr_C_0_2 corr_C_App if A==1 & uncensored_C==1,nocons
gen rca = _b[corr_C_App]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_App]				/*This is corr(u_C,u_A)*/
matrix semoments=semoments\_se[corr_C_App]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_C_0_TBUB		=_b[corr_C_0_2]
scalar rho_C_App_TBUB	=_b[corr_C_App]

biprobit (uncensored_C =$spec_U) (Aw =	$c_R),nolog	
predict a_unc3,xb1
predict a_Aw,xb2
scalar rho_unc_aw=e(rho)

egen c2=rmean(a_unc2 a_unc3)
gen  c3=A_App_3ml
gen  c4=A_Aw_3ml
scalar r23=rho_unc_ap
scalar r24=rho_unc_aw 
scalar r34=rho_aw_ap
scalar r34_2=(r34-r23*r24)/(sqrt(1-r23^2)*sqrt(1-r24^2))
scalar r24_3=(r24-r23*r34)/(sqrt(1-r23^2)*sqrt(1-r34^2))
scalar r23_4=(r23-r24*r34)/(sqrt(1-r24^2)*sqrt(1-r34^2))
matrix Sigma = (1, r23, r24 \ r23, 1, r34 \ r24, r34, 1)

mdraws, neq(3) dr(20) prefix(q) random
matrix C=cholesky(Sigma)
egen F3 = mvnp(c2 c3 c4), chol(C) dr(20) prefix(q) adoonly

gen corr_C_00=normalden(c2)*binormal((1-r23^2)^(-0.5)*(c3-r23*c2),(1-r24^2)^(-0.5)*(c4-r24*c2),r34_2)/F3
gen corr_C_ap=normalden(c3)*binormal((1-r23^2)^(-0.5)*(c2-r23*c3),(1-r34^2)^(-0.5)*(c4-r34*c3),r24_3)/F3
gen corr_C_aw=normalden(c4)*binormal((1-r24^2)^(-0.5)*(c2-r24*c4),(1-r34^2)^(-0.5)*(c3-r34*c4),r23_4)/F3

gen v=u-(rho_C_0_TBUB)*corr_C_00-(rho_C_App_TBUB)*corr_C_ap			/*use consistent estimates of corr.coeff. from other parts of the code*/

***TABLE 10: RHO_{C,R}
reg v corr_C_aw if Aw==1 & uncensored_C==1 & A==1,nocons
gen rcr = -_b[corr_C_aw]
scalar b_CA=_b[corr_C_aw]
scalar b_CAse=_se[corr_C_aw]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_aw]
matrix semoments=semoments\_se[corr_C_aw]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
eststo rc_coeff: estpost su rcl rca rcr,meanonly

esttab bigml1 bigml2 bigml3 using $figures\tableOA5.tex, ///
	cells(b(fmt(3)) se(par fmt(3)) mean(fmt(3))) title("ML Estimation") label unstack replace style(tex) longtable drop(c21 c31 c32 athrho cantwork_notemp: lnsigma) ///
	stats(N sample_avg_V sample_avg_A sample_avg_L sample_avg_R sample_avg_x het_V het_se_V het_L  het_se_L het_A het_se_A het_R het_se_R sigmac se_sigmac) 

esttab bigml1 bigml2 bigml3 using $figures\table6_1.tex, ///
	keep(female vign_f logbs) ///
	cells(b(fmt(3)) se(par fmt(3)) mean(fmt(3))) title("ML Estimation") label unstack replace style(tex) ///
	stats(N sample_avg_V sample_avg_A sample_avg_L sample_avg_R sample_avg_x het_V het_se_V het_L  het_se_L het_A het_se_A het_R het_se_R sigmac se_sigmac) 

	
matrix M=moments,semoments
matrix rownames M = 1_bF_V 2_bFv_V 3_corr_L_V 4_b_F_L 5_bF_A 6_corr_L_A 7_corr_R_A 8_corr_L_R 9_bF_C 10_sdu_C 11_EuC_L 12_EuC_A 13_EuC_Aw
matrix M=M[4..6,1..2]\M[1..3,1..2]\M[8..8,1..2]\M[9..11,1..2]\M[13..13,1..2]\M[12..12,1..2]\M[7..7,1..2]

preserve
	matrix list M
	svmat M
	keep M1 M2
	keep if M1!=.
	gen str20 namevar=""
	replace namevar="bF_L"            in 1
	replace namevar="bF_A"            in 2
	replace namevar="corr_L_A"        in 3
	replace namevar="bF_V" 			  in 4
	replace namevar="bFv_V"           in 5
	replace namevar="corr_L_V"        in 6
	replace namevar="corr_L_R"        in 7
	replace namevar="bF_C"            in 8
	replace namevar="sdu_C"           in 9
	replace namevar="EuC_L"           in 10
	replace namevar="EuC_Aw"          in 11
	replace namevar="EuC_A"   		  in 12
	replace namevar="corr_R_A"        in 13
	save $data\moments_for_R,replace
restore

clear

log close


